% this is the script for calculating the fractal dimention (D-value) in 2D
% Version : 1.2

clc;
clear all;
close all;

%Read circle data file
f_id=fopen('2D-section.dat');
XYR=fscanf(f_id,'%f %f %f\n',[3 inf]);
f_return=fclose(f_id);

%X=XYZR(1,:);
%Y=XYZR(2,:);
R=XYR(3,:);
D=R*2;

maxD=max(R)*2;
minD=min(R)*2;
%NumSphere=size(R,2);

%Calculate D-value
binInterval=0.1;
tempBin=[-10:binInterval:10];
maxDLog10=log10(maxD);
maxDLog10=min(find(tempBin>maxDLog10));
minDLog10=log10(minD);
minDLog10=max(find(tempBin<minDLog10));

phiBin=seqreverse([tempBin(minDLog10):binInterval:tempBin(maxDLog10)]);
numPhiBin=size(phiBin,2)-1;
Dref=1;
DBin=Dref*power(10,phiBin);

cumulCount=0;

for i=1:numPhiBin
    iBin=find(D<=DBin(i) & D>DBin(i+1));
    count=size(iBin,2);
    cumulCount=cumulCount+count;
    cumulFreq(i)=cumulCount;
    if count==0
        diameter(i)=(DBin(i)+DBin(i+1))/2;
    else
        diameter(i)=median(D(iBin));
    end;
end;

%Plot data
    hold on;

    logCFreq=log10(cumulFreq);
    logDiameter=log10(diameter);
    plot(logDiameter,logCFreq,'o');

    rangeInit=1;
    %rangeFinal=numPhiBin;
    rangeFinal=14;

    logCFreq=log10(cumulFreq(rangeInit:rangeFinal));
    logDiameter=log10(diameter(rangeInit:rangeFinal));
    
    plot(logDiameter,logCFreq,'o','MarkerFaceColor','b');
    %loglog(diameter(rangeInit:rangeFinal),cumulFreq(rangeInit:rangeFinal),'+');
    
    bestFit=polyfit(logDiameter,logCFreq,1);
    fitX=linspace(min(logDiameter),max(logDiameter),100);
    fitY=polyval(bestFit,fitX);
    plot(fitX,fitY,'--r','LineWidth',2);
    text(mean(fitX),mean(fitY),['y = ',num2str(bestFit(1)),'x','+',num2str(bestFit(2)),' \rightarrow'],'HorizontalAlignment','right','FontSize',13)
    
    axis_font_size=13;
    xlabel('Diameter (log)','fontsize',axis_font_size);
    ylabel('Cumulative Frequency (log)','fontsize',axis_font_size);
    
    hold off;

    D_Value=bestFit(1)*-1

